
################
## DID Method ##
################



colnames(wmat) = rownames(wmat)
wellmat = as.matrix(txPanel$cumuwells)
windmat = as.matrix(txPanel$cumucap)
urbanmat = as.matrix(txPanel$urban)
pair.shale <- matrix(nrow = 100,ncol = 5)
pair.wind <- matrix(nrow = 200,ncol = 5)
k = 1;
l = 1;
for (i in 1:254) {
  for (j in 1:254) {
    if (wmat[i,j] > 0 && any(wellmat[i,] > 10)&& all(wellmat[j,] == 0) && any(urbanmat[i,] == urbanmat[j,])) {
        pair.shale[k,] = c(index(wmat)[i],rownames(wmat)[i], index(wmat)[j], colnames(wmat)[j], urbanmat[j,1])
        k = k+1 
      }
    if (wmat[i,j] > 0 && any(windmat[i,] > 50)&& all(windmat[j,] == 0) && urbanmat[i,1] == urbanmat[j,1]) {
        pair.wind[l,] = c(index(wmat)[i],rownames(wmat)[i], index(wmat)[j], colnames(wmat)[j], urbanmat[j,1])
        l = l+1 
      }
  }
}

pair.shale = na.omit(pair.shale)
pair.wind = na.omit(pair.wind)

fm.shaledid <- emp ~ wells + shalecty + postshale + shalecty*postshale

coef.did <- matrix(nrow = dim(pair.shale)[1],ncol = 3)

k = 1;
for (i in 1:dim(pair.shale)[1]) {
txPair1 <- subset(txPanel, attr(txPanel,"index")[,1] == pair.shale[i,2])
txPair2 <- subset(txPanel, attr(txPanel,"index")[,1] == pair.shale[i,4])
txPair2[,"postshale"] = txPair1[,"postshale"]
txPair <- rbind(txPair1, txPair2)

emp.did = lm(fm.shaledid, txPair)

coef.did[k,1] = coef(emp.did)[4]
if (!is.na(coef.did[k,1])) {
coef.did[k,2] = summary(emp.did)$coefficients["shalecty:postshale","Pr(>|t|)"]
coef.did[k,3] = txPair[1,"urban"]
}
k = k+1
}
coefShale.Pvalue = coef.did[coef.did[,2]<0.001,]
coefShaleP = mean(coefS.Pvalue)
coefR.shale = coef.did[coef.did[,3] == 1,]
coefRS = mean(coefR.shale[,1])
coefU.shale = coef.did[coef.did[,3] == 2 & coef.did[,2]<0.001,]
coefUS = mean(coefU.shale[,1])

## wind
fm.winddid <- emp ~ newcap + windcty + postwind + windcty*postwind

coef.did <- matrix(nrow = dim(pair.wind)[1],ncol = 3)

k = 1;
for (i in 1:dim(pair.wind)[1]) {
txPair1 <- subset(txPanel, attr(txPanel,"index")[,1] == pair.wind[i,2])
txPair2 <- subset(txPanel, attr(txPanel,"index")[,1] == pair.wind[i,4])
txPair2[,"postwind"] = txPair1[,"postwind"]
txPair <- rbind(txPair1, txPair2)

## Graphs of pair.wind counties

## xyplot(emp~year, data = txPair.Wind, groups = attr(txPair,"index")[,2], type = "o", auto.key=list(title="County", space = "right", cex=1.0))

emp.did = lm(fm.winddid, txPair)

summary(emp.did)

coef.did[k,1] = coef(emp.did)[4]
if (!is.na(coef.did[k,1])) {
## coef.did[k,2] = summary(emp.did)$coefficients["windcty:postwind","Pr(>|t|)"]
coef.did[k,2] = summary(emp.did)$coefficients["windcty:postwind","Pr(>|t|)"]
coef.did[k,3] = txPair[1,"urban"]
}
k = k+1
}

coefW.Pvalue = coef.did[coef.did[,2]<0.001,]
coefWP = mean(coefW.Pvalue)
coefR.wind = coef.did[coef.did[,3] == 1,]
coefRW = mean(coefR.wind[,1])
coefU.wind = coef.did[coef.did[,3] == 2 & coef.did[,2]<0.001,]
coefUW = mean(coefU.wind[,1])
